      module lanczos_gamma
!     corrections in factorial calculation with double precision gamma built-in function
      use prec
      implicit none

      contains

    recursive function lancz_gamma_g7(a) result(g)
!   Gamma function by Lanczos method
    implicit none

    real(dp), intent(in) :: a
    real(dp) :: g
 
    real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp) 
    integer, parameter :: cg = 7
 
    ! these precomputed values are taken by the sample code in Wikipedia,
    ! and the sample itself takes them from the GNU Scientific Library
    real(dp), dimension(0:8), parameter :: p = &
         (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, &
         771.32342877765313, -176.61502916214059, 12.507343278686905, &
         -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)
 
    real(dp) :: t, w, x
    integer :: i
 
    x = a
 
    if ( x < 0.5_dp ) then
       g = pi / ( sin(pi*x) * lancz_gamma_g7(1.0_dp-x) )
    else
       x = x - 1.0_dp
       t = p(0)
       do i=1, cg+2
          t = t + p(i)/(x+dble(i))
       end do
       w = x + dble(cg) + 0.5_dp
       g = sqrt(2.0_dp*pi) * w**(x+0.5_dp) * exp(-w) * t
    end if
    end function lancz_gamma_g7

      function lancz_gamma_g9(z) result(reslt)
!     Lanczos approximation of Gamma function
!     Adopted from Matlab file of :
!     Paul Godfrey (pgodfrey@intersil.com)
!            References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
!            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
!            Y. Luke, "Algorithms ... functions", 1977
!            J. Spouge,  SIAM JNA 31, 1994. pp. 931
!            W. Press,  "Numerical Recipes"
!            S. Chang, "Computation of special functions", 1996
!
      implicit none
      real(dp) :: reslt
      real(dp), intent(in) :: z
      integer, parameter :: g = 9
      real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)
      real(dp), dimension(g+2), parameter :: c = &
      (/1.000000000000000174663, &
     5716.400188274341379136, &
   -14815.30426768413909044, &
    14291.49277657478554025, &
    -6348.160217641458813289, &
     1301.608286058321874105, &
     -108.1767053514369634679, &
        2.605696505611755827729, &
       -0.7423452510201416151527e-2, &
        0.5384136432509564062961e-7, &
       -0.4023533141268236372067e-8 /)
      integer :: k
      real(dp) :: s,ss,LogofGamma,t

      t=z+dble(g)
      s=0.0_dp
      do k=g+2,2,-1
        s=s+c(k)/t
        t=t-1.0_dp
      enddo
      s=s+c(1)
      ss=(z+dble(g)-0.5_dp)
      s=log(s*sqrt(2*pi)) + (z-0.5_dp)*log(ss)-ss

      LogofGamma = s
      reslt = exp(LogofGamma)

      end function lancz_gamma_g9

      end module lanczos_gamma
